package com.qp;

import com.joptimizer.functions.ConvexMultivariateRealFunction;
import com.joptimizer.functions.LinearMultivariateRealFunction;
import com.joptimizer.functions.PDQuadraticMultivariateRealFunction;
import com.joptimizer.optimizers.JOptimizer;
import com.joptimizer.optimizers.OptimizationRequest;
import org.apache.commons.math3.linear.Array2DRowRealMatrix;
import org.apache.commons.math3.linear.RealMatrix;

/**
 * Created with IntelliJ IDEA.
 * User: bigmar
 * Date: 12/11/14
 * Time: 2:22 PM
 * To change this template use File | Settings | File Templates.
 */
public class QPCaller
{
    public static double[][] call(double[][] cMatrix, double[] mixtureDistribution, double[] etaEstimate ) throws Exception
    {
        // Objective function
//        double[] qvector={2,-35, -47};

        double[][] cTransposeC=new double[cMatrix.length][cMatrix.length];
        for (int i=0; i<cMatrix.length; ++i)
            for (int j=0; j<cMatrix.length; ++j)
            {
                for (int k=0; k<cMatrix.length; ++k)
                    cTransposeC[i][j]+=cMatrix[k][i]*cMatrix[k][j];
            }
        double[] q=new double[cMatrix.length];

        // Build q vector
        for (int i=0; i<cMatrix.length; ++i)
        {
            for (int j=0; j<cMatrix.length; ++j)
                q[i]-=etaEstimate[j]*cMatrix[j][i];
        }


        PDQuadraticMultivariateRealFunction objectiveFunction = new PDQuadraticMultivariateRealFunction(cTransposeC, q, 0);

        int length=(int)Math.round(Math.sqrt(cMatrix.length));
        double[][] A = new double[length*2-1][cMatrix.length];
        double[] b = new double[length*2-1];

        // Build equalities sum(Aij)=1 for columns
        for (int i=0; i<length; ++i)
        {
            for (int j=0; j<length; ++j)
            {
                A[i][i+j*length]=1;
            }
            b[i]=1;
        }
        // Build stationary distribution equality
        for (int i=length; i<2*length-1; ++i)
        {
            int mixtureIndex=i-length;
            for (int j=0; j<length; ++j)
            {
                A[i][mixtureIndex*length+j]=mixtureDistribution[j];
            }
            b[i]=mixtureDistribution[mixtureIndex];
        }
        //inequalities
        ConvexMultivariateRealFunction[] inequalities = new ConvexMultivariateRealFunction[cMatrix.length];
        for (int i=0; i<length; ++i)
        {
            for (int j=0; j<length; ++j)
            {
                double[] parameters=new double[cMatrix.length];
                parameters[i*length+j]=-1;
                inequalities[i*length+j] = new LinearMultivariateRealFunction(parameters, 0);
            }
        }
//

        //optimization problem
        OptimizationRequest or = new OptimizationRequest();
        or.setF0(objectiveFunction);
        or.setFi(inequalities);
        or.setA(A);
        or.setB(b);
        or.setToleranceFeas(1.E-12);
        or.setTolerance(1.E-12);

        //optimization
        JOptimizer opt = new JOptimizer();
        opt.setOptimizationRequest(or);
        try
        {
        int returnCode = opt.optimize();
        }
        catch (Exception e)
        {
            System.out.println();
            return null;
        }
        double[] sol = opt.getOptimizationResponse().getSolution();

        double[][] retVal=new double[length][length];
        for (int i=0; i<length; ++i)
        {
            for (int j=0; j<length; ++j)
                retVal[i][j]=sol[i+j*length];
        }
        return retVal;
    }

    public static double[] call1(double[][] cMatrix) throws Exception
    {
//        for (int i=0; i<cMatrix.length; ++i)
//            for (int j=0; j<cMatrix.length; ++j)
//                if (cMatrix[i][j]<0)
//                    cMatrix[i][j]=-cMatrix[i][j];

        PDQuadraticMultivariateRealFunction objectiveFunction = new PDQuadraticMultivariateRealFunction(cMatrix, null, 0);

        //equalities
        double[][] A = new double[cMatrix.length][1];
        double[] b = {1};
//        double [] q=new  double[cMatrix.length];
        double [] initialPoint=new  double[cMatrix.length];
        for (int i=0; i<cMatrix.length;++i)
        {
            A[i][0]=1;
//            q[i]=2+Math.random();

            initialPoint[i]=Math.random();
        }

        //inequalities
//        ConvexMultivariateRealFunction[] inequalities = new ConvexMultivariateRealFunction[1];
//        inequalities[0] = new LinearMultivariateRealFunction(q, -1);

        //optimization problem
        OptimizationRequest or = new OptimizationRequest();
        or.setF0(objectiveFunction);
        or.setInitialPoint(initialPoint);
//        or.setFi(inequalities); //if you want x>0 and y>0
//        or.setA(A);
//        or.setB(b);
        or.setToleranceFeas(1.E-12);
        or.setTolerance(1.E-12);

        //optimization
        JOptimizer opt = new JOptimizer();
        opt.setOptimizationRequest(or);
        int returnCode = opt.optimize();
        System.out.println(returnCode);
        double[] sol = opt.getOptimizationResponse().getSolution();
        return sol;
    }
    public static void main(String[] args)
    {
        double[][] cMatrix=new double[2][2];
        cMatrix[0][0]=1;
        cMatrix[0][1]=2;
        cMatrix[1][0]=3;
        cMatrix[1][1]=4;


        RealMatrix matrix=new Array2DRowRealMatrix(cMatrix);
        RealMatrix transposeMatrix=matrix.transpose();
        RealMatrix multiplyMatrix=transposeMatrix.multiply(matrix);
        double[][] cTransposeC=new double[cMatrix.length][cMatrix.length];
        for (int i=0; i<cMatrix.length; ++i)
            for (int j=0; j<cMatrix.length; ++j)
            {
                for (int k=0; k<cMatrix.length; ++k)
                    cTransposeC[i][j]+=cMatrix[k][i]*cMatrix[k][j]/2;
            }
        System.out.println(multiplyMatrix);
        System.out.println(new Array2DRowRealMatrix(cTransposeC));
    }
}
